ExtractKBranchCombos <- function(BranchDepths,BranchProbs,CompProbs,K) {
  
  # This function is the intersection of N-Choose-K & Part's Hierarchical Rules.
  # It is designed to identify all the ways in which K can be summed from 
  # B Branches, each with different depths. 
  #
  # Inputs
  #   BranchDepths - C x 1 vector of integers
  #   BranchProbs  - C x 1 list vector, with each C having a 
  #                  vector length = BranchDepths[c] of probabilities for that event
  #   CompProbs - C x 1 vector of that component's overall reliability probability 
  #
  # Outputs
  #   ProbKFromBranch - scalar - 
  #
  # The algorithm is written as a state space tracker, where for K failures 
  # we are tracking each k (as in 1:K) as it's own "object" of sorts.  
  # The k-objects represent one count of the sub-component's branches.
  # These k-objects have rules about where they begin, how they interact 
  # with their neighbors, how they search their space, and when they can 
  # go no futher - algorithm has to end sometime! Those rules are:
  # K-OBJECT-CENTRIC
  #   1. All K of the k-objects begin immediately adjacent to their neighbors.
  #   2. The k-objects begin in the "left most" or first component branch.
  # COMPONENT CENTRIC
  #   3. For each Component Branch, the k-objects must start at the first position
  #   4. For each Component Branch, the k-objects must be adjacent to each other
  #   4.1 k-objects can only be disjoint or disconnected from each other by 
  #       residing in separate Component Branches
  # ALGORITHM CENTRIC
  #   5. The algorithm is done when all k-objects are adjacent at the last/end
  #      Component Branch
  #   6. The algorithm can have no more than (C^K)/2 + C^(K-1)/2 steps.
  #      This assumes that all Component Branch depths are greater than K,
  #      when this is not true, the hypercube analogy becomes a sparsely 
  #      populated cube.
  #
  # Written by: Sean M. Gonzalez, Willcor Inc., November 2012
  #=========================================================================
  
  # First, collect all the Component Branch information, i.e. the context, 
  # for all the k-objects
  numComp = length(BranchDepths)
  BranchStates = vector("list",numComp)
  for (c in 1:numComp) {
    BranchStates[[c]] = logical(BranchDepths[c])
    cat("Initialized BranchStates =",BranchStates[[c]],"\n")
  }
  
  # A useful value for position tracking of asymmetric component depths
  cumDepths = integer(numComp)
  for (c in 1:numComp) {
    cumDepths[c] = sum(BranchDepths[1:c])
  }
  if (any(BranchDepths==0)) { stop(BranchDepths) }
  
  #=========================================================================
  # "Stars and Bars" has a total of (C-1)choose(K-1) combinations
  permn(c(1,2,3,0,0))
  #=========================================================================
  # Second, create the k-objects. Each will have a set of states that need 
  # to be tracked after each action. Those states are:
  #   1.   Number
  #   2-3. Position - Component Branch index and depth in that branch
  #   4.   Step position (as opposed to branch position)
  #   5.   Aft k-object number (0 for the first k-object)
  #   6-7. Final position at algorithm end
  #   8.   Trigger increment of this k-object
  #   9.   total K
  kObjects = integer(K*9) # 9 states
  c(K,9) -> dim(kObjects) # I have K as the rows so I imagine them moving up 
                          #   and down the Component Branches
  for (k in 1:K) { # Set the states
    kObjects[k,1] = k
    # Next available Component Branch unoccupied position
    {
    kpos     <- mapStep2Branch(k,numComp,cumDepths)
    finalStep = cumDepths[numComp]-K+k
    finalPos <- mapStep2Branch(finalStep,numComp,cumDepths)
    cat("kpos =",kpos,"\n")
    BranchStates[[kpos[1]]][kpos[2]] = TRUE} # Now assign it
    
    kObjects[k,2] = kpos[1]
    kObjects[k,3] = kpos[2]
    kObjects[k,4] = k # The step
    kObjects[k,5] = k-1 # Will be zero for the first
    kObjects[k,6] = finalPos[1] # -max( min(0,(k+1)-cumDepths[numComp]) ,-(k+1))
    kObjects[k,7] = finalPos[2] # numComp # Do the assignment steps backwards!
    if (k != K) { kObjects[k,8] = FALSE } # Trigger
    else        { kObjects[k,8] = TRUE  }
    kObjects[k,9] = K
  }
  for (c in 1:length(BranchStates)) {
    cat("Populated BranchStates =",BranchStates[[c]],"\n")
  }
  cat("k-objects initialized! \n")
  #=========================================================================
  # Now that I have the states created and assigned/polulated, I have to
  # operate on those states, count the events, and accumulate the probabilities
  eventsCounted = 0
  eventProbs = 0
  allStatesEnded = FALSE
  while (!allStatesEnded) {
    
    # Assert (double check) that kObject positions add to K
    {
      testK = 0
      for (b in 1:numComp) {
        testK = testK + sum(BranchStates[[b]])
      }
      if (testK != K) {warning("Danger Will Robinson!! Danger!")}
      eventsCounted = eventsCounted +1 # One more in the Branched-N-Choose-K
    }
    
    # Now calculate the probability for these positions by:
    #   1. Multiplying P(Branch-Failures) *P(Branch-Survives)
    #   2. Adding each independent event probability
    {
      activeBranches = logical(numComp)
      xx = unique(kObjects[,2])
      activeBranches[xx] = TRUE
      thisProb = 1
      
      # Non-Active branches need to survive
      for (b in 1:length(BranchStates)) {
        cat("BranchStates =",BranchStates[[b]],"\n")
      }
      for (a in 1:numComp) {
        if (activeBranches[a]) {
          xx = sum(BranchStates[[a]])
          thisProb = thisProb * BranchProbs[[a]][xx] } # Prob X fail in branch
        else {
          thisProb = thisProb * (1-CompProbs[a])} # Prob whole branch fails
      }
    }
    # Events are independent and their probabilities add as such
    cat("thisProb =",thisProb,"\n")
    eventProbs = eventProbs + thisProb
    cat("eventProbs =",eventProbs,"\n")
    
    # I COULD CALCULATE THE COST OF THAT STATE, AND THIS IS THE LOWEST FIDELITY POSSIBLE
    # I COULD CALCULATE THE COST OF THAT STATE, AND THIS IS THE LOWEST FIDELITY POSSIBLE
    # I COULD CALCULATE THE COST OF THAT STATE, AND THIS IS THE LOWEST FIDELITY POSSIBLE
    
    # Finally, update all the k-objects, unless this is the last state.
    cat("Current State: ",kObjects[,2:3]," - End State:",kObjects[,6:7],"\n")
    allStatesEnded = all((kObjects[,2:3]-kObjects[,6:7])==0)
    if (allStatesEnded) {break}
    else {
      for (k in K:1) {
        # Need this object's next Position (nextPos)
        # The next position is determined by the 'odometer' approach to k-objects, 
        # where the 2nd k-object doesn't move until the 1st has gone through it's motions, 
        # 3rd until the 2nd, etc.
        updatedObjects <- stepKObj(k,K,kObjects[k-1,],kObjects[k,],BranchDepths,numComp,cumDepths)
        kObjects[k,]   = updatedObjects$thisObj
        kObjects[k-1,] = updatedObjects$aftObj
        
        # Update the branch states & k-object
        {
        Bx = kObjects[k,2]
        By = kObjects[k,3]
        BranchStates[[Bx]][By] = FALSE
        # Update the k-object
        kObjects[k,] = newObj
        Bx = kObjects[k,2]
        By = kObjects[k,3]
        BranchStates[[Bx]][By] = TRUE
        }
      }
#       stop()
    }
  }
  BNcK = list("NumEvents"=eventsCounted,"Prob"=eventProbs)
  return(BNcK)
}
#===============================================================================
stepKObj <- function(k,K,aftObj,thisObj,BranchDepths,numComp,cumDepths) {
  
  # This function is designed to apply the rules for incrementing each kObject
  # which ultimately results in the entire set of objects occupying the next event.
  # The rules for this process are:
  # 1. For any given Branch, all occupied position must be continuous and begin at 1
  # 2. The maximum positions of the kObjects in each branch must sum to K
  # 3. An object has gone through all it's potential states given it's aftObj 
  #      configuration when it either reaches its final position
  #      or is occupying it's final Branch and can't violate rule #1
  # 4. Each kObject's neighbor must be in the first or last position of the adjacent well
  
  # Only take steps if you're triggered, or you're the first object.
  if ((thisObj[1] != K) && (!thisObj[8])) {
    # Do nothing
    cat("DO NOTHING!?!?",thisObj,"\n")
  }
  # must comply with the 1-X (1 through X) rule
  # We have lowered the rope down the well as far as possible, now it must be the next well
  # IS THIS NECESSARY!!!!
  else if (thisObj[1]==K && K < BranchDepths[thisObj[2]]) {
    # if ((branch[1]+1) > numComp)
    # I don't need to account for the commented case above because 
    # I will update each of the k-objects in turn
    thisObj[2] = (thisObj[2] +1)
    if (thisObj[2] >numComp) {
      thisObj[2] = thisObj[2]%%numComp
    }
    thisObj[3] = 1
    # I do not need to increment the next object because 
    # I have not created a gap with that object.
  }
  # Since this function only looks at one object at a time
  # We have to check to see if this is the last object, 
  # and if it's in its final position.
  else if ((thisObj[1] == 1) && all(thisObj[6:7] == thisObj[2:3])) {
    # We're done!
    print("Last pony in the stable. Done!!")
    thisObj[8] = FALSE
  }
  # Apply Rule 3
  # Object is in its final position, or 
  # I've explored all the states given the aft object position
  else if (all(thisObj[6:7] == thisObj[2:3]) || 
          (thisObj[2]==thisObj[6] && (thisObj[4]-aftObj[4])>1 ) ) { 
    # Increment based on the aftObj next position and 
    # flag thisObj to trigger its aftObj
    cat("aftObj Before",aftObj[2:3]," -",numComp," -",cumDepths," - k=",k," - K=",K,"\n")
    aftObj[8] = TRUE
    thisObj[2:3] <- detNextStep(aftObj[2:3],numComp,cumDepths,2)
    cat("thisObj After",thisObj[2:3],"\n")
  }
  # Increment the object and
  else {
    # Increment this object and ensure K is the sum
    cat("thisObj Before -",thisObj[2:3]," -",numComp," -",cumDepths," - k=",k," - K=",K,"\n")
    thisObj[2:3] <- detNextStep(thisObj[2:3],numComp,cumDepths,1)
    cat("thisObj After",thisObj[2:3],"\n")
  }
  
  # Populate the output
  updatedObject = list("thisObj",thisObj)
  updatedObject$aftObj = aftObj
  return(thisObj)
}
#===============================================================================
detNextStep <- function(position,numComp,cumDepths,numSteps) {
   
#   nxt = list("Step"=NA,"Pos"=NA)
  nxt = integer(2)
#   kStep <- mapBranch2Step(position,cumDepths)
  if (position[1]==1) { cumPos = 0 }
  else {                cumPos = cumDepths[position[1]-1]}
  
  kStep = cumPos + position[2]
  nextStep = kStep +numSteps
  
  # Steps have to loop around, modular to branches
  cat("Next Step (before) = ",nextStep,"\n")
  if (nextStep > cumDepths[numComp]) {
    nextStep = nextStep%%cumDepths[numComp]
  }
  cat("Next Step (after) = ",nextStep,"\n")
  newPosition  = mapStep2Branch(nextStep,numComp,cumDepths)
  return(newPosition)
}
#===============================================================================
mapStep2Branch <- function(step,numComp,cumDepths) {
#   ativeBranch = logical(numComp)
#   cat("step=",as.numeric(step)," & cumDepths=",as.numeric(cumDepths),"\n")
  xx = which(step <= cumDepths)
#   activeBranch[xx[1]] = TRUE
  depth = max(1,step -cumDepths[xx[1]]) #cumDepths[activeBranch])
  branchPos = c(xx[1],depth)
  return(branchPos)
}
#===============================================================================
mapBranch2Step <- function(branch,cumDepths) {
  # branch - 1x2 vector giving position in the branches 
  #        - the output of mapStep2Branch
  step = cumDepths[branch[1]] + branch[2]
  return(step)
}